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' We perform Monte Carlo simulations on the hard-core attractive Yukawa system to test the 

04 . Optimized Baxter Model that was introduced in [P.Prinsen and T. Odijk, J. Chem. Phys. 121, 

}_( ■ 6525 (2004)] to study a fluid phase of spherical particles interacting through a short-range pair 

' potential. We compare the chemical potentials and pressures from the simulations with analytical 

predictions from the Optimized Baxter Model. We show that the model is accurate to within 10 
percent over a range of volume fractions from 0.1 to 0.4, interaction strengths up to three times 
' the thermal energy and interaction ranges from 6 to 20 % of the particle diameter, and performs 

CNJ , even better in most cases. We furthermore establish the consistency of the model by showing that 

the thermodynamic properties of the Yukawa fluid computed via simulations may be understood on 
the basis of one similarity variable, the stickiness parameter deflned within the Optimized Baxter 
I Model. Finally we show that the Optimized Baxter Model works significantly better than an often 

O , used, naive method determining the stickiness parameter by equating the respective second virial 

. ' coeflicients based on the attractive Yukawa and Baxter potentials. 
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I. INTRODUCTION 



o , 

I ^ I . In a recent paper [ll two of us devised a method to approximate systematically a system of spherical hard particles 
that interact through a short-range pair potential by a system of particles interacting via an effective Baxter potential 
. The latter consists of a hard-core repulsion and a sticky attraction at the surface of the particles which is computed 
^ ' by a variational principle for the free energy. The original potential was a sum of attractive and repulsive contributions 
\ (i.e. a square well plus a Debye-Hiickel interaction ^J) but the variational method also applies to a purely attractive 

I ■ interaction provided its range is sufficiently smaller than the particle diameter. The advantage of approximating the 

original interaction by the Baxter potential is that the fluid phase of the Baxter model has been studied extensively, 
both theoretically and in computer simulations 0, ^1 ■ This means that once the 

\^ correspondence between the two systems has been established, all the results of the Baxter model can be used for the 
: original system. 

In the Optimized Baxter Model (OBM) 0, the free energy of the actual system is functionally expanded in terms 
of the Mayer function where the reference state is a suspension of hard spheres interacting via an effective sticky 
y ■ potential. The stickiness parameter associated with the latter is determined by setting the first-order term in this 
' , \ expansion equal to zero. This constitutes a variational principle since the second-order term turns out to be either 
^ ■ positive or negative definite. Nevertheless, the exact nature of the expansion is difficult to assess analytically. For 
Q instance, there may be mathematical problems arising from the limiting procedure in which the range of the effective 
Q I adhesion goes to zero as its magnitude becomes infinitely large Thus, a computational test of the OBM is 



important. 

The model system we consider consists of hard-sphere particles with an attractive Yukawa interaction 



H ■ Uy{x) I oo 0<x<2 

^^1-/3-^^ ->2 



r 

X = —. 

a 

Here a is the radius of the particles, fcs is Boltzmann's constant, T is the temperature, (3 is the dimensionless well 
depth and a/n is a measure of the range of the attractive tail (if we set the actual well depth = unity, /3 may be 
viewed as identical with l/ZcsT). The somewhat awkward factors 1/2 in front of the dimensionless coordinate x are 
caused by the fact that we want to scale all distances by a as in Ref. |ll whereas, at the same time, we wish to use 
the same notation as in other papers on simulation studies of the same system |l4l Il5j where distances are scaled by 
the diameter 2a of the particles. 
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The liquid-solid coexistence of this system has been studied before at various values of k [ij, ll5| . These papers 
do not report the chemical potentials and pressures at coexistence however. We need these data to compare them 
to the predictions of the OBM to test its validity. We therefore perform new simulations to determine the volume 
fraction, chemical potential and pressure at various points along the phase boundaries. Moreover, we also determine 
the chemical potential and pressure within the fluid region of the phase diagram so as to gauge the accuracy of the 
OBM at lower concentrations.. 

We start by reviewing equations relevant to the OBM as applied to the Yukawa potential (Eq. JQ) in the next 
section. In section III we describe the numerical simulations which, in section IV, are compared with the theoretical 
predictions. 



II. THEORY 

Here we give a short summary of the theory developed in Ref. 1]. The relevant equations needed to determine 
the effective adhesion parameter t and some of the thermodynamic properties of the system are presented here. For 
details of the derivation we refer to Ref. and references mentioned there. 

We consider a system of spherical particles of radius a. The interaction U between the particles is pairwise additive 
and consists of a hard sphere repulsion plus a short-range interaction Ui that is either purely attractive or consists of 
a combination of attractive and repulsive interactions (range <^ a). In the latter case the attraction has to be strong 
enough to compensate for the repulsion — we come back to this issue later. For convenience, all distances are scaled 
by the radius a of the particles so we have 

. . ( CO 0< X <2 , , 

= x>2 (2) 

r 

X = —. 

a 

We wish to replace this system by a suspension of adhesive hard spheres with the same radius which is our reference 
state. The interaction of the latter is given by the adhesive hard sphere (AHS) potential of Baxter ^ 



Uahs{x) 



oo < X < 2 

Inf^ 2<x<2 + u:. (3) 
x>2 + u} 



Here, r is a positive constant that we wish to determine and which signifies the strength of the effective adhesion and 
the limit w | has to be taken after formal integrations. The reason for approximating the original system by the 
AHS system is that the latter has been conveniently solved in the Percus-Yevick approximation |2| • This means that 
once the correspondence between the two systems has been established by appropriately choosing t, other properties 
like for example the chemical potential, the pressure and the compressibility of the system can easily be computed 
analytically from the solution of the AHS system. 

We first describe how to choose the stickiness parameter r. In the limit of vanishing density this is done by equating 
second virial coefficients since we must equate the respective free energies of the two systems. 



This amounts to choosing 



Here, B^^ ~ IGTra'^/S is the second virial coefficient of a solution of hard spheres. At finite densities this procedure 
necessarily breaks down, however, because the higher virials come into play. The stickiness parameter r, which 
depends on the density, has to be obtained by identifying the free energy of the actual system with that of the 
reference state as well as possible. In the functional expansion of the excess free energy in terms of the Mayer function 
[T6| . we then demand that the first order correction vanishes. This leads to the condition 0] 

9 \ 

dxa;2(e-^(-)/'=-^-l)g(x) = ^, (6) 
3 
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where g{x) is the regular part of the pair correlation function g(x) of the reference AHS system and 

^ 6{r,+ il-rj)T) ( _ I rK2 + ^) \ 

,7(1 -r?) 1^ V 6{v+{l-v)rr) ^' 

with rj the volume fraction of particles. For x < 2, g{x) equals zero owing to the hard-core repulsion, whereas ^(x) 
tends to unity for large x. Since the interaction U (x) is of short range, we approximate g{x) in the interval 2 < a; < 4 
by the first two terms of its Taylor expansion 

X <2 

-2)) 2<a;<4. (8) 
X > 4 




Here, we define 



and 



G = At (9) 



77 fvi^-v).2 1 + 11??^ ,1 + 577 9(1 + 77) 1 



2r(l-77)V 12 12 1-77 2(1 -77)2 Ay ^ ' 

At a given volume fraction 77, t can then be determined iteratively from Eqs. l(H|l- H10() . An iteration scheme, which 
converges fast, consists of choosing a starting value of r, determining A from Eq. Q, then r from 

2A-3/rdxa;2(e-^(-)A«i^-l) ^ ^ 

r= 3 ^ , 11 

3Xj^dxx^il + H{x-2))ie-u(x)/ksT 

A again from Eq. Q and so on until convergence to the required accuracy is achieved. 

There are two cases in which the above method does not yield meaningful results. The first occurs when the 
short-range interaction has both attractive and repulsive components in the event that the repulsion is too strong 
in comparison with the attraction. The total interaction is then effectively repulsive in nature so it is clear that a 
suspension of particles interacting in such a way cannot be approximated by an AHS system. In this case, the iteration 
scheme described above leads to a r which keeps on increasing and does not converge. If tq is negative in the limit of 
vanishing density (Eq. ^) implying a net repulsion, it is advisable not to compute r also, even though it could lead 
to positive values at higher densities. Secondly, the attraction may be too strong. There exists a critical value of the 
stickiness parameter Tc below which there is a range of densities for which there is no real solution of A 

r,= '^. (12) 

This means that if the attraction is strong enough (i.e. when t is too small), there will not be a positive real solution 
to Eq. (|ll|l . In this case the iteration scheme would produce complex values of t. 

We now state several thermodynamic properties resulting from the solution of the Baxter model which we will 
compi: 

compressibility route 0, 13 

Pvo _ + V + V^) '7^(l + ^?/2), 

^ (1 - 77)3 (1-77)2 + 36 

and 

Here vq = Attu^ /3 is the volume of a particle, 



need further on. To compute the pressure P and the chemical potential we use the expressions derived via the 

A + ^A^ (13) 



A(l-77) 



6re(18rr,- 1)2 A(l - 77) - ISr^ 



l-6r. 



3 2 2 377(1 + 477) , , 677(2 + 77) I877 6(T-Te)2 
^-2'^^- (I-77) ^+Tr^"T^""r.(l-6r,)^" 

(15) 

is the contribution to the chemical potential that vanishes in the hard-sphere limit (r 00) and the chemical potential 
of the reference state (in the context of the Baxter model) is given by 



where h is Planck's constant and m is the mass of a particle. 
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III. SIMULATIONS 



We perform Monte Carfo simulations at constant volume V and temperature T on a system of = 256 hard 
spheres with a short-range Yukawa attraction so we have (compare with Eq. |Q) 



Ui{x) _ ^ e~«(^/2-i) 



x/2 



(17) 



We introduce a cutoff at a; = 4, so that U{x) = for x > 4. We determine the Helmholtz free energy per particle 
/at at a chosen set of parameters of /3, k and rj by thermodynamic integration at constant k and rj, starting from the 
known free energy per particle of the hard sphere system (/3 — 0) which is defined at the same volume fraction (see 
e.g. [11) 



fNiv,P)-fNiv,0) 



(18) 



N 



Here {Ui)n is the average energy per particle where the average is computed in the state with f3 = /3' . From this we 
determine the equation of state z^irj, 13) 



pksT 



pksT 



d_ 

^ drj 



fNiv,f3)-fN{v,0) 



kuT 



(19) 



where the particle density p is related to the volume fraction by pvQ — rj, and the chemical potential pn(j]jP) is given 
by 



Pn{v^P) Mjv(?7,0) 9 



knT 



knT 



drj 



fN{v,/3)-fNiv,0) 
knT 



(20) 



The expression for the pressure is exact for a system consisting of a finite number of particles N, whereas that for 
the chemical potential has an error of order N~'^, because in the simulations we change only the volume V of the 
box leaving the number of particles invariant. (See the appendix for details). For the equation of state of the pure 
hard-sphere system we have 



z'(r,,0) = 1 + 



47] + 1.2162247]^ + 1.246720r/3 



1 - 2.19594477-1- 1.210352ry2 
valid when the system is a fluid 0| . It is quadrature to determine the chemical potential 

p\f^,0) 



knT 



= lm]-l + z''(t],0) + / 

^0 



For the pressure of the hard sphere (fee) solid we use [T3 

3 0.5921(^77-0.7072) 



z^(77,0) 



and for the chemical potential 

M^(^?,0) 
knT 



In 77 -1 + z%ri,0) + 



rim) 

knT 



",^,.^(V,0)-1 



(21) 



(22) 



(23) 



(24) 



where f'^{rji^)/kBT — 5.91889(4) is the free energy of the hard sphere solid in the thermodynamic limit TV —> 00 at 
volume fraction 7/0 — 0.5450 [23 • We thus calculate the pressure and the chemical potential of the system in the 
thermodynamic limit using 



^^^Vz (^,0)+77- 



Ui 
knT 



N 



and 



ksT keT drj^Ja ' \A:bT/ 



(25) 



(26) 



TV 
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where i = s,l. These expressions are not exact but correct to order because the number of particles in the 
simulations is finite (see the appendix for details). 

To determine the average energy per particle (C/i/fc^T) n we need to initiate the simulation by choosing a convenient 
starting configuration. In the case that the system is a solid, we assume it is an fee crystal at the appropriate density. 
For the liquid, a configuration at the required density is initiated by putting the particles in the box at random and 
then running the simulation until the particles no longer overlap. This is done at a low value of f3 and we then use 
this starting configuration for all values of /? at the same density. The simulation is then run for 10,000 cycles (i.e. 
trial moves per particle) at the relevant value of /3 to determine the appropriate maximum displacement of a particle 
at an acceptance probability of a particle displacement of 0.40. The maximum displacement is then fixed and the 
simulation is run for another 10,000 cycles for the system to equilibrate. Finally, the average energy per particle is 
measured every 100 cycles during another 50,000 cycles. 

To perform the integration in Eq. (|18|l we run simulations at values of /3 ranging from 0.1 until the appropriate 
value at intervals of 0.1. A simulation at /3 = 0.02 is also performed. We then fit the points to a curve and use this 
to perform the integration. To determine the density dependence of the free energy /jv about a certain density, we 
compute the free energy at about 10 values of the density close to it, at intervals of 0.1. We again fit these to a curve 
which is used in Eqs. and to determine the chemical potential and the pressure at the desired density. 

IV. RESULTS AND DISCUSSION 
A. Phase equilibrium 

We first test the Optimized Baxter Model on the fluid-solid coexistence of hard spheres with Yukawa attraction. 
This coexistence has been studied before via computer simulations |l4lll5j| but these papers did not report the pressure 
and chemical potential, data we do need here. 

At a given strength (3 and inverse range n of the attraction, we compute the volume fractions of the coexisting fluid 
and solid phases by equating the pressures and the chemical potentials in the respective phases. This is done for k = 7 
and 9, and /3 running from to 2 at intervals of 0.25 (see Table I). Our phase boundaries at k = 7 agree well with 
those computed by Dijkstra j3| where the same method was used though with a smaller system of = 108 particles. 
The deviation in volume fraction is at most 2% (we determined the phase boundaries from a plot presented in Ref. 
[THj l , so this may account for part of the difference) . At low values of /3 the agreement with the simulations of Hagen 
and Frenkel [lj| is also good, but with increasing /3 the difference between their phase boundary on the fluid side and 
ours becomes appreciable until our prediction of the volume fraction is about 20% higher than theirs at /3 = 2. We 
note that in Ref. 14] a different method was used to determine the phase boundary. The phase boundaries on the 
solid side do agree within 3%. We regain essentially the same picture at k = 9 though the difference in the phase 
boundaries at the fluid side is less pronounced (about 14% at /? = 2). The phase diagram at k = 9 was not determined 
in Ref. 

Next, we use the OEM to determine the effective stickiness parameter r (Eq. (|ll|l ) and the coexistence curve. By 
way of comparison, we also evaluate tq by equating the respective second virial coefficients of the attractive Yukawa 
interaction and the Baxter potential (see Eq. (|SJ)) and computing the properties of the resulting Baxter fluid. We 
will refer to this as the B2 method which is strictly correct only at very low concentrations. We employ Eqs. H13() 
and H14|) to calculate the pressure and the chemical potential from the volume fractions and the respective values of 
T from the two methods. These predictions are compared with the simulations in Fig. 1 (k = 7) and Fig. 2 (k = 9). 
It is clear from the figures that the predictions of the OBM are significantly better than those via the B2 method 
along the whole phase boundary. The OBM is actually quite accurate to within a few percent. Recall that at (3 — 0, 
i.e. in the absence of attraction, the two volume fractions predicted by the two methods necessarily coincide simply 
because t — 00 m that case. However, they do not agree with the simulations which is due to the fact that we use 
the accurate equation of state (Eq. 1211) ') in the simulations whereas the analytical theory is of course approximate. 
The latter overestimates the pressure and the chemical potential in the case of hard spheres. 

B. Consistency test 

We next check the consistency of the OBM. The Baxter model itself has been solved in the Percus-Yevick ap- 
proximation and we here use the compressibility route to obtain the thermodynamic properties. We know, however, 
that in the case of the hard sphere system, the analytical calculations via the Percus-Yevick approximation and the 
compressibility route are too high (e.g. at 77 = 0.4, both the pressure and the chemical potential are overestimated by 
4%). We therefore seek to test the argumentation within the OBM in a way which is less sensitive to approximations 
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inherent in the Baxter modeL For instance, we note that the stickiness parameter r in the OBM merely depends on 
the properties of the distribution function g very close to the sphere (see Eq. ^). Though this does depend on the 
Percus-Yevick approximation, it stands to reason that the functions G and H are more robust than the oscillatory 
behavior which Tj actually displays in full (and which is implicit in Eqs. and (|14|l 'l. Thus, in the following 

simulations, we investigate whether similarity is achieved with respect to the parameter r as given by Eq. 

Our procedure is as follows. We start at a given volume fraction. Next we choose a set of values of the inverse range 
of the Yukawa potential (i.e. k = 5, 7, 9, 11, 13 and 15). We then fix a certain value of the stickiness parameter t and 
compute the concomitant value /3 for each k with the help of Eq. Hll|) . If similarity does apply, the thermodynamic 
properties should depend solely on r and ry i.e. they ought to be independent of n at constant r. 

We have performed this test on simulations in a suitable range of volume fractions 77 and stickiness parameters r 
with associated interaction parameters k and f3 as chosen above. (See Figs. 3 to 6.) In some cases the attraction 
is so strong in terms of /3 that the simulated fluid is actually in a metastable region. In effect, if the system were 
macroscopic, phase separation into fluid and crystal phases would occur. (See for example the relevant figures in Refs. 
[T^ This happens when r = 0.1, for example in the case n = 15 where we have P ~ 3.5 for all four volume 

fractions and in the case k — 7 where f3 ~ 2.4. Despite the pre-emption of phase separation, we may still determine 
the pressure and chemical potential as if the phases were stable. 

We first note that the simulated thermodynamic properties are generally quite independent of k. (See the filled 
symbols in Figs. 3 to 6). This implies that r is indeed a useful similarity variable and the OBM is a consistent 
approximation scheme. The variation in the pressures and chemical potentials computed by simulation is only a few 
percent with few exceptions (e.g. in some instances, at lower values of k when the volume fraction is 0.4, or when 
r = 0.1 at the higher densities; in the latter case, scrutiny of the simulation snapshots shows that gelation seems to 
be occurring — note that the attraction is so strong that we are now well beyond the percolation threshold 12]). 

Next, it is of interest to compare the magnitudes of the simulated thermodynamic properties with those computed 
with the help of the OBM (see the curves in Figs. 3 to 6 which are horizontal because r was forced to be constant in 
each case). The analytical predictions are virtually quantitative, except at those densities at r — 0.1 where gelation 
seems to occur as discussed above and with regard to some of the pressures at higher concentrations. The latter are 
overestimated at r = 0.5 and 1 in Figs. 5 and 6 which we attribute to deficiencies in the Baxter model itself, since 
the simulational data are quite independent of k as stressed above. 

For the sake of comparison we have also displayed thermodynamic properties computed by the B2 method. At a 
certain n and /3 we evaluate tq with the help of Eq. |SJ| using the Yukawa interaction Eq. Q (thus tq is not constant 
like t) and then calculate the pressure and chemical potential within the Baxter model. The B2 method works well 
at low concentrations (see Fig. 3), which is not surprising since neglecting to variationally adjust virials higher than 
second is not so crucial in this case. However, the B2 method worsens progressively as the concentration increases 
and ultimately becomes unreliable (see Figs. 4-6). This is of course expected: the B2 method merely adjusts a single 
coefficient B2 whereas the free energy itself is variationally optimized in the OBM. 

We conclude that the Optimized Baxter Model is a convenient quantitative, analytical theory for computing the 
thermodynamic properties of a fluid of hard spheres interacting by an attraction of short range. Moreover, the 
variational scheme used in deriving the OBM is consistent, especially when the range of the potential is short i.e. less 
than approximately 10 % of the particle diameter {k > 10). Overall, the OBM is accurate to within 10 percent, except 
under conditions of very strong attraction at high volume fractions (r — 0.1, r/ = 0.3 and 0.4), and it is actually much 
more precise in most cases. 
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V. APPENDIX A 

Here we show that the error incurred in Eq. H20|) for the chemical potential of the iV-particle system is of order 
A^~^, whereas Eq. (|19|l for the equation of state is exact. We also prove that the error in the free energy of the system 
is of order N~^. 

Our simulations are carried out at a constant number of particles N. Hence, we modify the volume fraction r/ by 
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altering the volume of the simulation box. The free energy difference per particle 



(27) 



is determined as a function of the volume fraction, so in effect it is a function of 77 (or p) and N (and of course (3 and 
At). The exact equation of state zn^tj, (3) for the A^-particle system is then 



pksT 



dV 



= zn{v,0) +7]- — - — , 



where F^^q, (3) — NfN{ri, /?) and the exact chemical potential is 



dN 



V,T 



dr] 



kuT 



1 d AfN{r,,P) 
N dN-^ knT 



(28) 



(29) 



Here, and in the rest of the appendix, we have switched to the new independent variables rj and N so that derivatives 
with respect to r/ are taken at constant N and derivatives with respect to N are taken at constant 77. We see from 
Eqs. and (^5)1 that Eq. has an error of order whereas Eq. ((T^ is exact. 

We now assume that we may Taylor expand A/Ar(r;,/3) for small values of at constant volume fraction. It's 
not obvious that this is allowed. In the case of a crystal for example, the first-order correction to the free energy 
per particle due to the fact that the number of particles is finite, is of order iV~^lnA^ [20I I2H. This correction is 
the same for systems of identical numbers of particles however, regardless of the interaction. Since our is the 
difference in the free energies per particle pertaining to the two respective crystals (with different pair potentials), the 
0{N~^ IniV) corrections simply cancel. Moreover, from Ref. we know that the leading higher order corrections 
to the free energy per particles are of order iV~^. These deliberations are confirmed in Figs. 7 and 8 which show 
that the leading corrections to the average dimensionless energy per particle (t/i/fcsr)jv are indeed of order N^^ 
at the representative values /3 = 1, k = 15 and p{2a)^ = 0.4 {rj = n/15 0.20944) for the fluid and p{2a)^ = 1.2 
[f] — 7r/5 « 0.62832) for the solid. Therefore, we conclude that the free energy per particle in a system containing an 
infinite number of particles is given by 



In the same manner, the equation of state Zoo{rj^[3) is then 



Zooiv, P) = Zo.{71, 0) + 7,- + O 



and the chemical potential is expressed by 



ksT ksT dri ksT \N 



(30) 



(31) 



(32) 



VI. APPENDIX B 



We estimate the second-order correction to the free energy (see Appendix C of Ref. 1] ) 



4 ' 



(33) 



This correction leads for instance to a correction to the dimensionless pressure Pvo/ksT approximately equal to 
— 277A Ij. The first part of the analysis in Appendix C of Ref. ,lj is also useful here and we again approximate Y by 



where 



r ~ ^ ( 9G + lOGH - 12 + ^ 



B{x) = g{x) 



dt tB{t) 



/ dt tB{t) I ds s^B{s) 

J2 \ \_J2 



exp 



Uy{x) Uahs{x) 



keT 



keT 



- 1 



(34) 



(35) 
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and A, G and H are given by Eqs. (f7|)- (|10|l . We split the pair distribution function g{x) in the reference state into 
guj(,x) and a regular part 'g{x) given by Eq. © (see also Ref. 0): 



g{x) = g{x) + g^{x) (36) 



with 



x<2 
= ^ ^^ + 0(1) 2<x<2 + uj (37) 
x>2 + bj. 

We then insert the expressions for the potentials Eqs. Q and ^ into Eq. H35|l and derive in the limit — ^ 

Ax xB{x) = -- + G I dx x{l + H{x - 2)){e-^^^''^'^^^ - 1) + 0(6--^) (38) 

and 



2 



/ dxx^B{x)=-— + G da;a;3(i + ^(2;_2))(e-^^(")/'=«^- l) + 0(e-'=). (39) 

In both cases the integration on the right hand side should run from a; = 2 to a; = 4, so extending the integrals to 00 
only introduces errors of order e^'^. In the OEM, r is determined by the condition that the first-order correction to 
the free energy vanishes 

/ dxx^B{x) = -— + G da;a;2(l + iJ(a;-2))(e-^^(")/'=^^-l) + O(e-«) = 0. (40) 
J2 "J J2 

This expression is used to rewrite Eqs. H38() and (|39|l 

) 1 poo 

dx xB{x) = -7:G dx x{x - 2)(1 + H{x - 2))(e~^^(^)/'=^^ - 1) + 0(e~'=), (41) 

^ Jo 



dx x'^Bix) =G j dx a;2(x - 2)(1 + H{x - 2))(e-^^(")/^-«^ - 1) + ©(e"") (42) 

which are readily approximated. We substitute y = cxp[— K(a;/2 — 1)] which ultimately leads to 

f°° AG 
/ dxxB{x)^—Ji{P)^0{K-'^) 
J2 



(43) 



and 



Here we have introduced 



/ da;x3B(x) =.-^Ji(/3)+0(At-3). (44) 

J2 



Ji(/3) = - / dy Iny. (45) 

An approximation for Ji(/3) that is accurate to within 1.4% in the relevant range < /3 < 3.52 is given by 

1 iR\ f P+y"^ 0</3<0.8 

^i(/^)- j 2.8l(eO-34/5-l) 0.8 </3< 3.52. ^^^^ 

Finally, we insert Eqs. (|43|l and H44|) into Eq. H34|) . We thus obtain an approximation for the second-order correction 
to the free energy 

A ^ (9G + lOGH - 18 + -j JUPW . (47) 
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In Table II we present typical values of A. The corrections to the pressure are very small (compare with Figs. 1 and 
2). 
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TABLE I: Volume fraction of particles in the coexisting fluid and solid phases as a function of /3 and k determined 
by the simulations. The stickiness parameter is computed via the Optimized Baxter Model (OBM) and the B2 method 
{B2). 
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TABLE II: Estimates of the second-order correction A to the free energy for rj — 0.4 at various r and k. A is 
determined from Eq. (|47|l within the approximation given by Eq. H46|l . 



Figure Captions and Figures 

FIG. 1. Dimensionless chemical potential and dimensionless pressure as a function of the strength /? of the Yukawa 
potential for the coexisting fluid and solid phases. Here k = 7. The diamonds and the fitted line are results from the 
simulations. The squares are predictions from the Optimized Baxter Model (at the same densities) and the triangles 
have been computed by the B2 method. 

FIG. 2. Same as Fig. 1 but now at k = 9. 

FIG. 3. Dimensionless chemical potential and dimensionless pressure as a function of the inverse range k of the 
Yukawa potential at volume fraction rj — 0.1. The solid symbols are results from the simulations, the horizontal lines 
are predictions from the Optimized Baxter Model at a variety of fixed values of r. In the simulations the strength 
of the attraction f3 is chosen in such a way that the Optimized Baxter Model gives the appropriate value of r: grey 
filled diamonds t — 1, grey filled squares r = 0.5, black filled triangles r = 0.2, black filled squares r — 0.15 and black 
filled diamonds r = 0.1. The corresponding open symbols have been computed by the B2 method. 

FIG. 4. Same as Fig. 3 but now at volume fraction 77 = 0.2. 

FIG. 5. Same as Fig. 3 but now at volume fraction 77 ~ 0.3. 

FIG. 6. Same as Fig. 3 but now at volume fraction 77 = 0.4. 

FIG. 7. Example of the dependence of the average dimensionless energy per particle {Ui/kBT) n in the fluid on the 
size of the system. Here k = 9, /3 = 1 and 77 — 7r/15 « 0.20944 {p{2aY = 0.4). N denotes the number of particles. 
FIG. 8. Same as Fig. 7 but now for the solid at 77 = 7r/5 « 0.62832 {p{2af = 1.2). 
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Figure 6 
































